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In this communication, we propose a metliod for obtaining isolated excited states within the 
Full Configuration Interaction Quantum Monte Carlo framework. This method allows for stable 
sampling with respect to collapse to lower energy states and requires no uncontrolled approxima- 
tions. In contrast with most previous methods to extract excited state information from Quantum 
Monte Carlo methods, this results from a modification to the underlying propagator, and does not 
require explicit orthogonalization, analytic continuation, transient estimators or restriction of the 
Hilbert space via a trial wavefunction. Furthermore, we show that the propagator can directly yield 
frequency-domain correlation functions and spectral functions such as the density of states which 
are difficult to obtain within a traditional Quantum Monte Carlo framework. We demonstrate this 
approach with pilot applications to the neon atom and beryllium dimer. 



Almost all of the many variants of projector Quan- 
tum Monte Carlo (QMC) rely on the properties of the 
operator e~^^ , where due to its similarity to the time- 
evolution operator, the variable /? is denoted 'imaginary 
time'. Generally, this imaginary time is discretized, and 
the operator iteratively applied as a short-time propaga- 
tor, in order to simulate its action in the large /3 limit 
Expressing an initial wavefunction in the eigenbasis of 
the Hamiltonian of interest, the application of this e~^^ 
propagator results in a projection onto the i^^ eigenstate 
proportional to e~^^' , where Ei is the energy of this 
eigenstate. It is clear to see that in this large /3 limit, 
the projection onto the lowest energy eigenvector domi- 
nates the wavefunction, whereby ground state properties 
can be extracted, assuming some overlap with the ini- 
tial wavefunction. While this formalism is clearly pow- 
erful, by construction it exponentially quickly projects 
out excited states of the system which may be of inter- 
est, and are of critical importance in the simulation of 
finite-temperature properties, reaction dynamics, photo- 
chemistry and many other areas. 

To date, isolating excited states of systems via projec- 
tor QMC methods has only been practical with a restric- 
tion on the projection to sample a space which is approxi- 
mately orthogonal to those of the lower energy states, via 
nodal constraint or orthogonalization against them 
in a subspace projection method^, More indirectly, 
statistical methods have been used on short periods of 
imaginary-time in order to isolate individual decay rates 
in the spectrum by analytic continuation to a real-time 
dynamicjil, Q. However, these approaches are not en- 
tirely satisfactory; accurate nodal surfaces for excited 
states can be difficult to obtain, resulting in a larger fixed 
node error and potentially transient estimators, while the 
subspace projection method has limited applicability 
In addition, the Bayesian techniques which rely on max- 
imizing entropy to approximate a notoriously unstable 
inverse Laplace transform, have difficulty achieving quan- 
titative accuracy within noisy data sets Despite 



this, there are examples of accurate excited states within 
the nodal constraint [llj , while maximum entropy tech- 
niques are particularly prevalent in solid state calcula- 
tions to obtain the density of states, often in the case of 
continuous-time QMC as applied to quantum impurity 
models and dynamical mean- field theory [l2, 13 1. 

Here, we take a different approach to the calculation of 
excited states, within the context of Full Configuration 
Interaction Quantum Monte Carlo (FCIQMC)0-ll6|. 
This recently introduced method applies the imaginary- 
time evolution propagator to a stochastic 'walker' rep- 
resentation of the wavefunction expressed in the full 
space of Slater determinants constructed from a single- 
particle basis of size M. Although this reintroduces 
a basis set error compared to those methods operat- 
ing in real space, it confers various advantages which 
mitigate this. The discrete basis allows for an efficient 
walker annihilation algorithm, which can exactly over- 
come the fermion sign problem in the sampling, pro- 
vided enough walkers are usedjlTf. The 'initiator' ap- 
proximation was formulated to maintain a high annihi- 
lation rate, and control the growth of noise in a system- 



atically improvable fashion|15l. Il6j. This has allowed far 



larger systems to be treated at an accuracy compara- 
ble to that of Full Configuration Interaction (FCI or ex- 
act diagonalization) , within small and systematically im- 
provable error bars. In addition, a semi-stochastic adap- 
tation of the algorithm [ISLas well the introduction of 
a partial nodal constraint [l9j and ideas from quantum 
chemistrv[2ol| hold promise of increased accuracy and ef- 
ficiency of the method. 

In order to project out a targeted excited state rather 
than the ground state, we propose the use of a projection 
operator of the form 



P{H) = e 



(1) 



For sufficiently large /3, this Gaussian propagator will 
result in the dominant eigenstate being the one closest 
in energy to the chosen value of the diagonal offset S, 
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termed the shift. In the eigenbasis of H, {\^i),Ei} with 
^'o representing the ground state, and starting from an 
initial wavefunction \iPt) with S ~ E^., it can be seen 
that the long time propagation results in 



cx lim 



E 



-l3\Ei-Ekf 



I*, 



(2) 

We note here that a projector of this form was proposed 
back in 1983 within continuum QMC approaches j2l[ [2^ . 
although due to sign problems, and significantly larger 
timestep errors resulting from the fact that is more 
singular than H, only one-electron systems were re- 
ported, and no modern implementation exists in the lit- 
erature. This issue of the timestep highlights another 
advantage of working in a finite basis, in that the spec- 
trum is bounded both from below and above. This allows 
for linearization of the short-time propagator, 
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without becoming unbound and dominated by very high 
energy states oscillating in time, and without incur- 
ring timestep errors in the final wavefunction so long 
as the timestep is less than an upper bound given by 
T < -g '~E~' Repeated application of the short- 
time propagator therefore results in a 'power-method' for 
states on the interior of the spectrum, rather than at the 
extrema, with similarities to filter diagonalization jisili^ . 

Propagation with Eq. ([T]) leads to a theoretical decay 
of state j from i as ((^i-S) -(Ej-S) )^ contrast 
with the ground state propagator, this rate of decay de- 
pends on S, with it being advantageous to choose S to 
be as close to the energy of the state of interest as pos- 
sible. However, even if S is chosen exactly, the projec- 
tion of the non-dominant states is slower compared to 
the ground state propagation, and we will return to this 
issue later. In addition, unless S is chosen exactly, the 
long-time propagation of the dynamic will result in a con- 
tinued projection onto a decaying function of all states, 
including the dominant one. For this reason, the fac- 
tor of A is introduced into the short-time propagator, 
such that at convergence, its value can be varied in order 
to maintain a constant normalization of the domi- 
nant wavefunction and a stable number of walkers. This 
is analogous to the variation of S in the ground state 
projection[14|. 

The differential formulation of the exact dynamic gov- 
erned by this propagator for a given component of the 
wavefunction, Ci, can be written as 

da 



dT^ 



{A-l)a~eAj2iH,j- 



-6ijS){Hjk~SjkS)Ck, (5) 



where e — as A — 1, and the application of has 
been decomposed by a resolution of the identity over the 
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FIG. 1: Convergence of the propagation to the second ex- 
cited state of He2 at 2.5 A separation in a cc-pVDZ basis. 
S was fixed at -3.65Eh, and A at 1.004 until 10,000 walkers 
were present, denoted by the dotted line, where A was varied 
in order to keep this number constant. After variation, the 
average value oi A — 1 was 1.4(6) xlO~^. 



connecting space of determinants j. This formulation is 
now amenable to stochastic integration with a discrete 
walker representation of the determinantal wavefunction 
coefficients C. As with the ground state projection, there 
is no unique stochastic algorithm for this dynamic, but 
the one which we found to be most efficient involves a 
double spawning cycle, which requires little additional 
overhead compared to the ground state algorithm, and 
no additional memory requirements. Each iteration, the 
determinants represented by k are run through, and a 
spawning step attempted to determinant j, in the same 
fashion as the ground state propagation, but in nega- 
tive time. This results in a spawning probability to a 
connected determinant j with a stochastically realised 
signed amplitude of -p^^, where P{k\j) represents the 
normalised probability to randomly select symmetry- 



connected determinant j from i. Successfully spawned 
walkers are subsequently propagated again in the same 
iteration with a now forwards-time signed amplitude of 
— -p^jjjy. Care must be taken that for determinant fc, the 
diagonal 'death' processes from the first application of H 
are now interpreted as spawning events, which are also 
subsequently propagated via Hij. Each iteration, the 
factor of A is applied initially as a separate enhancement 
of the local population of each determinant, with the 
absolute population on each determinant growing with 
probability ACk- 

In Fig. [TJ we present an illustrative example of the al- 
gorithm for the helium dimcr in a cc-pVDZ basis, small 
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FIG. 2: Convergence of the projected energy estimate to 
the exact eigenvalue. Tlie reference determinant for the pro- 
jection was dynamically adjusted to project onto the largest 
weighted determinant in the sample. 



enough such that the full spectrum of cigenstates can be 
calculated and the convergence of the method analysed. 
The value of S was fixed at ^ 40mEh higher than the 
second excited state, but such that it remained the dom- 
inant state in the dynamic, and was subsequently found 
to be correctly projected out over time. This is despite 
working in a canonical representation, and starting with 
a single walker on the Hartree-Fock determinant which 
had an initial overlap with the ground state of close to 
one. In order to grow the walkers, A was initially fixed 
at a value of 1.004, and was varied when a target number 
of walkers was reached, in common with the procedure 
for the ground state propagation. The convergence of the 
projected energy, as defined in Ref. |14| . is shown in Fig. [5] 
for the same simulation, and reflects the decay from the 
sampled wavefunction of the first excited state. 

In order to reliably extend to larger molecular systems, 
it is worth considering how to transfer the salient ele- 
ments of the initiator approximation into this new prop- 
agation. The basis of this approximation is to attempt 
to propagate walkers corresponding to wavefunction sig- 
nal exactly, while walkers judged to be potentially noise 
are propagated with a truncated Hamiltonian which acts 
only over the instantaneously occupied subspace[l5l.[l6|. 
This is systematically improvable as the instantaneously 
occupied subspace grows, or the criterion for walkers cor- 
responding to signal becomes more inclusive. The sepa- 
ration between walkers corresponding to signal and po- 
tential noise is not unique. However, it seems sensible to 
retain the tested feature from ground state propagation 
that newly spawned walkers on previously unoccupied 



determinants (i) at the end of an iteration, must have 
come from an initial determinant {k) which is deemed 
to have a well-established sign, and therefore a popula- 
tion of walkers above nadd- Consequently, all walkers 
from the application of the first Hamiltonian operator 
are kept, while the information as to whether Ci is above 
the initiator criterion is passed through to the annihila- 
tion stage of the final set of spawned walkers. No walkers 
are therefore aborted over the resolution of the identity 
between (determinants j in Eq. (O). 

To test this on a larger system, we study an ex- 
cited state deep in the spectrum of the 10-electron neon 
atom in a cc-pVDZ basis, with an energy of approx- 
imately 2.5Eh above the ground state. Setting S to 
equal the CISDTQ energy for the corresponding state, 
and while remaining in a canonical Hartree-Fock ba- 
sis and starting from a random distribution of walkers 
throughout the whole space, we achieved a converged 
energy of -126.2118(4)Eh, compared to the FCI value 
of -126.21177Eh. This value is 4.76mEh lower than the 
initial guess provided by CISDTQ. It would be highly 
advantageous to develop a robust algorithm for varying 
S dynamically during the run, as is done for the ground 
state algorithm. This could be used to maximise the rate 
of growth of walkers or alternatively minimize A, both of 
which should adjust S to more closely match the eigen- 
value of the state, remove reliance on the initial guess 
and increase the convergence rate. However, since this 
requires finding a minimum in a quadratically varying 
and noisy dataset, no robust algorithm has been identi- 
fied so far. 

Dynamic correlation and response functions due to 
some perturbation, either in the frequency or time do- 
main, are of critical importance in electronic structure 
theory [2^, and are directly measured in experiments 
to probe the electronic properties of materials through 
techniques such as neutron scattering or photoelectron 
spectroscopy [2^. Many methods, including in general 
QMC approaches, have significant difficulty in calculat- 
ing these quantities ^3] 5 often having to rely on unsta- 
ble analytic continuation from imaginary time correla- 
tion functions [sl. [gI. [sl-fioj. while other methods can bias 
towards high or low energy regimes [l^. Although other 
correlation functions are possible, here we look at the ex- 
ample of an advanced Green's function, a central concept 
in electronic structure where the 'perturbation' at time 
i = is the creation of a hole in orbital j. For nega- 
tive time periods, i, these can be written in the time and 
frequency domain respectively as 
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uj-{H-Eo) + iS 



(6) 
(7) 



Unlike the inverse Laplace transform required for the an- 
alytic continuation of imaginary time correlation func- 
tions, the transform between these two domains is a well- 
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conditioned and numerically stable fourier transform in 
the presence of noisy data. Spectral density functions, 
such as the density of states for extended systems, are 
then defined in the Lehmann representation [28| as 



1 



which in the small S limit tends to 



(8) 
^(9) 
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where 6 in the above equation represents the dirac-dclta 
function. 

Assuming A = 1, application of the propagator in 
Eq. ([T]) for a time = ^yr will result in the wavefunction 



(11) 



which when applied to an initial wavefunction |?/'t) = 
fljl^'o) obtained from the ground-state dynamic, and 
then projected onto --^(^'Q|aJ will result in the distri- 
bution 
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for = w + Eq . This will tend to the spectral function 
given in Eq. pUj) in the large limit. The real parts of 
the Green's function can then by obtained if needed from 
the Kramers-Kronig relation[25[. We note that a related 
Green's function can be obtained directly by integrating 
Eq. over /3, with the addition of the small imaginary 
component S to the dynamic. Unfortunately however, 
this integral is only convergent for {lu — E^~^ + Eq) > S, 
and so the FCIQMC calculation will blow up at the poles. 
In systems with a continuous spectra, this would not be 
appropriate, and so we do not pursue this approach here. 
Results from a pilot investigation of the beryllium dimer 
in a cc-pVDZ basis, where the exact Green's function can 
be obtained from complete diagonalization, are shown in 
Fig.El 

In order to reduce the statistical error, it may be nec- 
essary to average over a small number of independent 
calculations at each frequency, and this can be combined 
with an elimination of the bias derived from choosing a 
correlated sample of (vPolal and ajlvPo) 01, by taking the 
^0 samples on each side of Eq. ([T2|) from different snap- 
shots in imaginary time. In addition, by storing multiple 
wavefunctions of the type (vPolal at the same time, all 
AP single-particle Green's functions can be calculated 




FIG. 3: High energy window of the spectral function 
A~{l,l,Lj) for exact propagation with 5 = 0.0141Eh, and 
stochastic evaluation via FCIQMC for an equivalent time 
P = 50a.u. for frozen-core Be2 in a cc-pVDZ basis at 2.5A . 
Vertical lines indicate the difference between the ground state 
energy and the eigenvalues of the N-1 system symmetry con- 
nected in G~(l, although some are coupled too weakly 
to contribute significantly to the spectral function. Approx- 
imately 10 independent calculations at each value of u; were 
averaged to obtain the errorbars. 



at a cost of ©[A/] FCIQMC calculations per frequency 
point, without the expectation of any variation in accu- 
racy between high and low energy regimes. 

However, despite modest successes, it is clear that ob- 
taining converged results through the use of this operator 
is substantially more difficult than with the ground state 
projection. This is mainly due to an additional factor of 
(tAE)~^ in the number of iterations required to project 
out undesired states with energy gap AE for comparable 
accuracy to the ground state propagation. The result is 
that while in the ground state propagation excited states 
were filtered relatively quickly with only isolated conver- 
gence issues in the case of near degeneracy [l^, the num- 
ber of iterations required for excited state propagation 
are substantially increased, as well as the dynamic being 
less well-conditioned with respect to walker fluctuations. 
This is also exacerbated by a generally more multicon- 
figurational wavefunction which increases random error 
in the projected energy estimator 1J|. A more judicious 
choice of orbital basis and initial conditions optimized 
for the state of interest, as well as a multireference pro- 
jected energy formulation(l8| would ameliorate many of 
these issues. In addition, there is the possibility of pre- 
conditioning techniques familiar from iterative diagonal- 
ization methods [29} being transferred into the stochastic 
dynamic, as well other operators, such as e~'^'^', which 
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should behave more efficiently and allow for extension to 
larger systems. Research in these directions is currently 
under way. It is clear that alternative propagators within 
the FCIQMC dynamic holds promise for obtaining accu- 
rate excited states. 
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